In [1]:
import sys
sys.path.insert(0, '/home/swl/Owncloud/Code/Edinburgh/Multiple Variables/')
# sys.path.insert(0, '/home/swk/Owncloud/Code/Edinburgh/Multiple Variables/')
In [2]:
''' Imports first etc. '''
import holoviews as hv
%load_ext holoviews.ipython
%output widgets='embed' holomap='widgets' size=200 fig='png'
import pylab as pl
import saveload as sl
import maxL as ml
from matplotlib.colors import LinearSegmentedColormap
from   pylab import exp,cos,sin,pi,tan
/usr/local/lib/python2.7/dist-packages/IPython/nbformat.py:13: ShimWarning: The `IPython.nbformat` package has been deprecated. You should import from nbformat instead.
  "You should import from nbformat instead.", ShimWarning)
/usr/local/lib/python2.7/dist-packages/IPython/nbconvert.py:13: ShimWarning: The `IPython.nbconvert` package has been deprecated. You should import from ipython_nbconvert instead.
  "You should import from ipython_nbconvert instead.", ShimWarning)
/usr/local/lib/python2.7/dist-packages/IPython/kernel/__init__.py:13: ShimWarning: The `IPython.kernel` package has been deprecated. You should import from ipykernel or jupyter_client instead.
  "You should import from ipykernel or jupyter_client instead.", ShimWarning)
In [3]:
## Holoviews settings and options
# set labels
# angles = pl.linspace(-90,90,65)
# rads   = 
# labels =[(i,int(angles[i])) for i in pl.linspace(0,64,3)]

# dimensions
p_dim = hv.Dimension('Preferred orientation',unit='deg')
a1_dim = hv.Dimension('Stim 1',unit='deg')
a2_dim = hv.Dimension('Stim 2',unit='deg')
r_dim = hv.Dimension('Response',unit='sp/s')
L_dim = hv.Dimension('Log Likelihood')
bias_dim = hv.Dimension('Bias',unit='deg')
In [4]:
%opts Overlay [show_frame=False show_grid=False]
%opts Curve [show_grid=False]
%opts Scatter [show_grid=False show_frame=False]
%opts Raster [colorbar=True]
# yticks=labels xticks=labels] (cmap='rainbow')

Coding and decoding of multiple variables in the brain

ANC Seminar

Sander Keemink

2015-09-22

In [5]:
reload(ml)
# function for simple firing rate
def fun(x,var=0,par = [1,1]):
    ''' 2 d tuning curve circular variables
    A neurons tuning is simply given by multiplication of two von mises functions
    
    Parameters
    --------------
        x: array
            x values (preferred orientations)
        var: double
            orientation
        par : double
            Parameters. Given as [a1,k1], the inputs into von mises function.
            
    Returns
    ---------------
    Array of firing rates given the var and par values
    '''
    return ml.mises(x,var,par[0],par[1])
# set bounds for fitting
bnds = [(-pi*1.5,pi*1.5)]

# set starting conditions for fitting
var0 = [-pi/2,pi/2]

# number of values/neurons
N = 64

# number of center orientations
Nc = 16

# model vars
kc = 3
ac  = 10

# basic x values for circular variables
A = pl.linspace(-pi,pi,N+1) #  in radiants
Ax= A*180/pi # in degrees for plotting purposes

# centre orientation values for holomap
C = pl.linspace(-pi,pi,Nc+1)

# create empty holomaps
hmap_f = hv.HoloMap(kdims=['Centre Orientation']) # underlying firing rate
hmap_fest = hv.HoloMap(kdims=['Centre Orientation']) # underlying firing rate
hmap_r = hv.HoloMap(kdims=['Centre Orientation']) # noisy response
hmap_vl = hv.HoloMap(kdims=['Centre Orientation']) # vertical lines

# fill holomaps
for c in C: # for a
    # get underlying response
    f = fun(A,c,par = [ac,kc])
    
    # get poisson response
    r = pl.poisson(f)
    
    # get decoded variables through likelihood maximization
    L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par=[ac,kc],noise='poisson')

    # get estimated underlying firing rate
    f_est = fun(A,var=var_max,par=[ac,kc])
    
    # add to holomap
    hmap_f[c*180/pi] = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'underlying response')
    hmap_fest[c*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response')
    hmap_r[c*180/pi] = hv.Scatter(zip(Ax,r),vdims=[r_dim],kdims=[p_dim],label = 'measured response')
    hmap_vl[c*180/pi] = hv.VLine(c*180/pi)
    

Coding

  • Simple tuning curve model
$$ f(\theta_c) = a\exp(k\cos[\phi-\theta_c]) $$
In [6]:
%%opts Scatter (alpha=0) Curve (alpha=0.5) Overlay [xticks=[-180,0,180] yticks=[0,10] show_legend=False]
showfig = hmap_f*hmap_r
In [7]:
showfig
Out[7]:
In [8]:
%%opts Scatter (alpha=1) Curve (alpha=0.5) Overlay [show_legend=True]
showfig = hmap_f*hmap_r
  • Measured firing rate $$ r = \texttt{Poisson}(f) $$
In [9]:
showfig
Out[9]:
  • What was the encoded orientation?

Maximum Likelihood Decoding

  • Likelihood to measure a firing rate $$ p(\mathbf{r}|c) = \prod_i p(r_i|c) $$

  • How to decode? asymptotically optimal thing is max likelihood decoding * $$ \hat{\theta}_c = \texttt{max}_{\theta_c} p(\mathbf{r}|\theta_c) $$

*Kay (1998) - Fundamentals of Statistical Signal Processing

In [10]:
# create empty holomaps
hmap_f = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Estimated underlying firing rate
hmap_fest = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Estimated underlying firing rate
hmap_r = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Estimated underlying firing rate
hmap_L = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Likelihood
hmap_Lloc = hv.HoloMap(kdims=['Estimated Centre Orientation']) # Likelihood

# set c
c = 0

# get underlying response
f = fun(A,c,par = [ac,kc])
curve_f = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'underlying response')

# get poisson response
r = pl.poisson(f)
scatter_r = hv.Scatter(zip(Ax,r),vdims=[r_dim],kdims=[p_dim],label = 'measured response')

# initiate empty likelihood array
L = pl.zeros(len(C))

# get Log Likelihood across estimated centre orientations
for i,c in enumerate(C): # for a    
    # get estimated underlying firing rate
    f_est = fun(A,var=c,par=[ac,kc])
    
    # get likelihood
    L[i] = ml.findLL(f,fun,A,c,par=[ac,kc])
    
    # add to holomap
    hmap_f[c*180/pi]    = curve_f
    hmap_r[c*180/pi]    = scatter_r
    hmap_fest[c*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response')
    hmap_Lloc[c*180/pi] = hv.Scatter(zip([c*180/pi],[L[i]]),vdims=[r_dim],kdims=[p_dim])(style={'s':100,'alpha':0.5})
    

curve_L = hv.Curve(zip(C*180/pi,L),vdims=[L_dim],kdims=[p_dim])

for i, c in enumerate(C):
    hmap_L[c*180/pi] = curve_L
In [11]:
%%opts Scatter (alpha=1) Curve (alpha=0.5) Overlay [xticks=[-180,0,180] yticks=3 show_legend=True]
showfig = hmap_f*hmap_r*hmap_fest+hmap_L*hmap_Lloc*hmap_Lloc
In [12]:
showfig
Out[12]:
In [13]:
%%opts Scatter (alpha=1) Curve (alpha=0.5) Overlay [xticks=[-180,0,180] yticks=[0,10] show_legend=True]
showfig = hmap_f*hmap_fest*hmap_r
In [14]:
showfig
Out[14]:
  • 'Optimal' variance (cramer-rao bound)
  • Non biased
  • I will be focussing on the bias

But neurons rarely respond to only single stimulus

Other examples:

  • overlapping stimuli
  • other variables/contexts (color, motion, speed, etc.)

Fig: Schwartz, Sejnowski & Dayan (2009)

Adjust coding model

$$f(\theta_1,\theta_2) = a\exp(k_1\cos[\phi-\theta_1])\{1-b\exp(k_2\cos[\phi-\theta_2])\} $$
In [15]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
    ''' 2 d tuning curve circular variables
    A neurons tuning is simply given by multiplication of two von mises functions
    
    Parameters
    --------------
        x: array
            x values (preferred orientations)
        var: array
            value of the two inputs
        par : array
            Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
            
    Returns
    ---------------
    Array of firing rates given the var and par values
    '''
    return ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(x,var[1],par[2],par[3]))

# set parameters
par  = [10,3,0.5,2]
par1 = [10,3,0,1] # effect of second stimulus turned off
par2 = [1,0,0.5,1] # effect of first stimulus turned off

# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))

# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]

# number of neurons
N = 32

# number of contexts
nA = 16

# iterations 
iterations = 500

# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))

# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)

# stim 1
a1 = 0

# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)

# predefine holomaps
hmap_scatters    = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1   = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2   = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f           = hv.HoloMap(kdims=['context']) # population response
hmap_g           = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h           = hv.HoloMap(kdims=['context']) # stim 2 modulation

# loop over stim 2's
for i,a2 in enumerate(A2):
    print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
    # set variables
    var=[a1,a2]
    
    # get underlying firing rate
    f = fun(A,var=var,par=par)
    hmap_f[a2*180/pi] = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'response')
    
    # get drive and modulation
    g = fun(A,var=var,par=par1)
    hmap_g[a2*180/pi] = hv.Curve(zip(Ax,g),vdims=[r_dim],kdims=[p_dim],label = 'Stim 1 drive')
    h = fun(A,var=var,par=par2)
    hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
    
    # do noise realizations
    print 'Realization (max: ' + str(iterations) + '):'
    for it in range(iterations):
        print '\r'+str(it),
        # get noisy firing rate
        r = pl.poisson(f)
        
        # do the maximization
        L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
    
        # store variables
        a1_est[i,it] = var_max[0]
        a2_est[i,it] = var_max[1]
    
    # make sure everything is within -pi,pi
    a1_est[i,a1_est[i,:]<-pi]+=2*pi
    a1_est[i,a1_est[i,:]>pi] -=2*pi
    a2_est[i,a2_est[i,:]<-pi]+=2*pi
    a2_est[i,a2_est[i,:]>pi] -=2*pi
    
    # get bias
    a1_bias[i] = ml.circstats(a1_est[i,:])[0]
    a2_meanest = ml.circstats(a2_est[i,:])[0]
    
    # compare estimate for a2 with actual a2 to make circular
    a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
    a2_bias[a2_bias<-pi]+=2*pi
    a2_bias[a2_bias>pi] -=2*pi
    
    # fill in relevant holomaps
    hmap_scatters[a2*180/pi]    = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
    hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    
# fill in illusion curves
for i,a2 in enumerate(A2):
    hmap_illusion1[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
    hmap_illusion2[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16
Realization (max: 500):
499 Doing stimulus 1 of 16
Realization (max: 500):
499 Doing stimulus 2 of 16
Realization (max: 500):
499 Doing stimulus 3 of 16
Realization (max: 500):
499 Doing stimulus 4 of 16
Realization (max: 500):
499 Doing stimulus 5 of 16
Realization (max: 500):
499 Doing stimulus 6 of 16
Realization (max: 500):
499 Doing stimulus 7 of 16
Realization (max: 500):
499 Doing stimulus 8 of 16
Realization (max: 500):
499 Doing stimulus 9 of 16
Realization (max: 500):
499 Doing stimulus 10 of 16
Realization (max: 500):
499 Doing stimulus 11 of 16
Realization (max: 500):
499 Doing stimulus 12 of 16
Realization (max: 500):
499 Doing stimulus 13 of 16
Realization (max: 500):
499 Doing stimulus 14 of 16
Realization (max: 500):
499 Doing stimulus 15 of 16
Realization (max: 500):
499
In [16]:
showfig = hmap_f*hmap_g*hmap_h
In [17]:
showfig
Out[17]:
In [18]:
%%opts Curve {+axiswise} Scatter {+axiswise}
showfig  = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g*hmap_h
In [19]:
showfig.cols(2)
Out[19]:
In [20]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig  = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g*hmap_h
showfig += hmap_scatters
In [21]:
showfig.cols(2)
Out[21]:
In [22]:
labels =[(i,int(A[i]*180/pi)) for i in pl.linspace(0,len(A)-1,3)]
/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:1: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  if __name__ == '__main__':
In [23]:
%%opts Raster [colorbar=True yticks=labels xticks=labels] (cmap='rainbow')
%%opts Scatter (alpha=0.25 color='k')
# neruongs
nA = 32
A  = pl.arange(-pi,pi,2*pi/nA)

# sampling
nS = 32
As = pl.arange(-pi,pi,2*pi/nS)

# set variables
varid = 11
var = [0,A2[varid]] #a1,a2

# get underlying firing rate
f = fun(A,var=var,par=par)
curve_f = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'underlying response')

# predefine holomaps
hmap_fest1= hv.HoloMap(kdims=['stim2']) # underlying population response
hmap_fest2= hv.HoloMap(kdims=['stim2']) # underlying population response
hmap_f   = hv.HoloMap(kdims=['stim2'])    # estimated population response
hmap_pos1 = hv.HoloMap(kdims=['stim2'])  # position of stimulus
hmap_pos2 = hv.HoloMap(kdims=['stim2'])  # position of stimulus

hmap_L   = hv.HoloMap(kdims=['stim2'])  # likelihood map

# get likelihood across estimations
L = pl.zeros((len(A),len(A)))
for i,a1 in enumerate(A):
    for j,a2 in enumerate(A):
        # get likelihood
        L[i,j] = ml.findLL(f,fun,A,[a1,a2],par)
    
a1a = 0
a1b = -pi/16
for i,a2 in enumerate(As):
    hmap_f[a2*180/pi]   = curve_f
    # get underlying firing rate situation 1
    f_est    = fun(A,var=[a1a,a2],par=par)
    hmap_fest1[a2*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response 1')
    # situation 2
    f_est    = fun(A,var=[a1b,a2],par=par)
    hmap_fest2[a2*180/pi] = hv.Curve(zip(Ax,f_est),vdims=[r_dim],kdims=[p_dim],label = 'estimated response 2')
        


# get normalized version of relevant scatter
var1_norm = len(A)*(a1_est[varid,:]+pi)/(2*pi)
var2_norm = len(A)*(a2_est[varid,:]+pi)/(2*pi)
scatter = hv.Scatter(zip(var2_norm,var1_norm))

showfig = hmap_f*hmap_fest1*hmap_fest2

Ambiguity in encoding model

In [24]:
showfig
Out[24]:
In [25]:
showfig = hv.Raster(L,kdims=[a2_dim,a1_dim])*scatter+hv.Empty()
In [26]:
showfig
Out[26]:

Stimulus-1-dependence

  • from
$$f_a(\theta_1,\theta_2) = a\exp(k_1\cos[\phi-\theta_1])\{1-b\exp(k_2\cos[\phi-\theta_2])\} $$
  • to
$$f_b(\theta_1,\theta_2) = a\exp(k_1\cos[\phi-\theta_1])\{1-b\exp(k_2\cos[\theta1-\theta_2])\} $$
In [27]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
    ''' 2 d tuning curve circular variables
    A neurons tuning is simply given by multiplication of two von mises functions
    
    Parameters
    --------------
        x: array
            x values (preferred orientations)
        var: array
            value of the two inputs
        par : array
            Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
            
    Returns
    ---------------
    Array of firing rates given the var and par values
    '''
    return ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(var[0],var[1],par[2],par[3]))
# set parameters
par  = [10,3,0.5,2]
par1 = [10,3,0,1] # effect of second stimulus turned off
par2 = [1,0,0.5,1] # effect of first stimulus turned off

# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))

# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]

# number of neurons
N = 32

# number of contexts
nA = 16

# iterations 
iterations = 500

# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))

# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)

# stim 1
a1 = 0

# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)

# predefine holomaps
hmap_scatters    = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1   = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2   = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f           = hv.HoloMap(kdims=['context']) # population response
hmap_g           = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h           = hv.HoloMap(kdims=['context']) # stim 2 modulation

# loop over stim 2's
for i,a2 in enumerate(A2):
    print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
    # set variables
    var=[a1,a2]
    
    # get underlying firing rate
    f = fun(A,var=var,par=par)
    hmap_f[a2*180/pi] = hv.Curve(zip(Ax,f),vdims=[r_dim],kdims=[p_dim],label = 'response')
    
    # get drive and modulation
    g = fun(A,var=var,par=par1)
    hmap_g[a2*180/pi] = hv.Curve(zip(Ax,g),vdims=[r_dim],kdims=[p_dim],label = 'Stim 1 drive')
    h = fun(A,var=var,par=par2)
    hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
    
    # do iterations
    print 'iteration (max: ' + str(iterations) + '):'
    for it in range(iterations):
        print '\r'+str(it),
        # get noisy firing rate
        r = pl.poisson(f)
        
        # do the maximization
        L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
    
        # store variables
        a1_est[i,it] = var_max[0]
        a2_est[i,it] = var_max[1]
        
        # by convenction set sign properly (since can't be told by decoder)
        # which 
        if pl.sign(a2_est[i,it]) != pl.sign(a2):
            a2_est[i,it] *= -1
    
    # make sure everything is within -pi,pi
    a1_est[i,a1_est[i,:]<-pi]+=2*pi
    a1_est[i,a1_est[i,:]>pi] -=2*pi
    a2_est[i,a2_est[i,:]<-pi]+=2*pi
    a2_est[i,a2_est[i,:]>pi] -=2*pi
    
    # get bias
    a1_bias[i] = ml.circstats(a1_est[i,:])[0]
    a2_meanest = ml.circstats(a2_est[i,:])[0]
    
    # compare estimate for a2 with actual a2 to make circular
    a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
    a2_bias[a2_bias<-pi]+=2*pi
    a2_bias[a2_bias>pi] -=2*pi
    
    # fill in relevant holomaps
    hmap_scatters[a2*180/pi]    = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
    hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    
# fill in illusion curves
for i,a2 in enumerate(A2):
    hmap_illusion1[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
    hmap_illusion2[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16
iteration (max: 500):
499 Doing stimulus 1 of 16
iteration (max: 500):
499 Doing stimulus 2 of 16
iteration (max: 500):
499 Doing stimulus 3 of 16
iteration (max: 500):
499 Doing stimulus 4 of 16
iteration (max: 500):
499 Doing stimulus 5 of 16
iteration (max: 500):
499 Doing stimulus 6 of 16
iteration (max: 500):
499 Doing stimulus 7 of 16
iteration (max: 500):
499 Doing stimulus 8 of 16
iteration (max: 500):
499 Doing stimulus 9 of 16
iteration (max: 500):
499 Doing stimulus 10 of 16
iteration (max: 500):
499 Doing stimulus 11 of 16
iteration (max: 500):
499 Doing stimulus 12 of 16
iteration (max: 500):
499 Doing stimulus 13 of 16
iteration (max: 500):
499 Doing stimulus 14 of 16
iteration (max: 500):
499 Doing stimulus 15 of 16
iteration (max: 500):
499
In [28]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig  = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g*hmap_h
showfig += hmap_scatters
In [29]:
showfig.cols(2)
Out[29]:

Multiple populations

$$f_a(\theta_1,\theta_2) = a_a\exp(k_{a1}\cos[\phi-\theta_1])\{1-b\exp(k_{a2}\cos[\phi-\theta_2])\} $$$$f_b(\theta_1,\theta_2) = a_b\exp(k_{b1}\cos[\phi-\theta_2])\{1-b\exp(k_{b2}\cos[\phi-\theta_1])\} $$
In [30]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
    ''' 2 d tuning curve circular variables
    A neurons tuning is simply given by multiplication of two von mises functions
    
    Parameters
    --------------
        x: array
            x values (preferred orientations)
        var: array
            value of the two inputs
        par : array
            Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
            
    Returns
    ---------------
    Array of firing rates given the var and par values
    '''
    f = pl.zeros(2*len(x))
    f[:len(x)] = ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(x,var[1],par[2],par[3]))
    f[len(x):] = ml.mises(x,var[1],par[0],par[1])*(1-ml.mises(x,var[0],par[2],par[3]))
    return f

# set parameters
par  = [10,3,0.5,2,10,3,0.5,2]

# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))

# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]

# number of neurons
N = 32

# basic x values for circular variables
A = pl.linspace(-pi,pi,N+1) #  in radiants
Ax = pl.linspace(-pi,3*pi,2*N+1)
Ax= Ax*180/pi # in degrees for plotting purposes

# number of contexts
nA = 16

# iterations 
iterations = 500

# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))

# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)

# stim 1
a1 = 0

# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)

# predefine holomaps
hmap_scatters    = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1   = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2   = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f           = hv.HoloMap(kdims=['context']) # population response
hmap_g           = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h           = hv.HoloMap(kdims=['context']) # stim 2 modulation

# loop over stim 2's
for i,a2 in enumerate(A2):
    print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
    # set variables
    var=[a1,a2]
    
    # get underlying firing rate
    f = fun(A,var=var,par=par)
    hmap_f[a2*180/pi] = hv.Curve(zip(Ax[:N],f[:N]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 1')
    
    # get second pop
#     g = fun(A,var=var,par=par1)
    hmap_g[a2*180/pi] = hv.Curve(zip(Ax[N:],f[N:]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 2')
    
    # get modulation
    h = fun(A,var=var,par=par2)
    hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
    
    # do iterations
    print 'iteration (max: ' + str(iterations) + '):'
    for it in range(iterations):
        print '\r'+str(it),
        # get noisy firing rate
        r = pl.poisson(f)
        
        # do the maximization
        L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
    
        # store variables
        a1_est[i,it] = var_max[0]
        a2_est[i,it] = var_max[1]
        
        # by convenction set sign properly (since can't be told by decoder)
        # which 
#         if pl.sign(a2_est[i,it]) != pl.sign(a2):
#             a2_est[i,it] *= -1
    
    # make sure everything is within -pi,pi
    a1_est[i,a1_est[i,:]<-pi]+=2*pi
    a1_est[i,a1_est[i,:]>pi] -=2*pi
    a2_est[i,a2_est[i,:]<-pi]+=2*pi
    a2_est[i,a2_est[i,:]>pi] -=2*pi
    
    # get bias
    a1_bias[i] = ml.circstats(a1_est[i,:])[0]
    a2_meanest = ml.circstats(a2_est[i,:])[0]
    
    # compare estimate for a2 with actual a2 to make circular
    a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
    a2_bias[a2_bias<-pi]+=2*pi
    a2_bias[a2_bias>pi] -=2*pi
    
    # fill in relevant holomaps
    hmap_scatters[a2*180/pi]    = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
    hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    
# fill in illusion curves
for i,a2 in enumerate(A2):
    hmap_illusion1[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
    hmap_illusion2[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16
iteration (max: 500):
499 Doing stimulus 1 of 16
iteration (max: 500):
499 Doing stimulus 2 of 16
iteration (max: 500):
499 Doing stimulus 3 of 16
iteration (max: 500):
499 Doing stimulus 4 of 16
iteration (max: 500):
499 Doing stimulus 5 of 16
iteration (max: 500):
499 Doing stimulus 6 of 16
iteration (max: 500):
499 Doing stimulus 7 of 16
iteration (max: 500):
499 Doing stimulus 8 of 16
iteration (max: 500):
499 Doing stimulus 9 of 16
iteration (max: 500):
499 Doing stimulus 10 of 16
iteration (max: 500):
499 Doing stimulus 11 of 16
iteration (max: 500):
499 Doing stimulus 12 of 16
iteration (max: 500):
499 Doing stimulus 13 of 16
iteration (max: 500):
499 Doing stimulus 14 of 16
iteration (max: 500):
499 Doing stimulus 15 of 16
iteration (max: 500):
499
In [31]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig  = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g#*hmap_h
showfig += hmap_scatters
In [32]:
showfig.cols(2)
Out[32]:

Asymetric populations

In [33]:
### context modulation, fixed
def fun(x,var=pl.array([0,0]),par=pl.array([1,1,1,1])):
    ''' 2 d tuning curve circular variables
    A neurons tuning is simply given by multiplication of two von mises functions
    
    Parameters
    --------------
        x: array
            x values (preferred orientations)
        var: array
            value of the two inputs
        par : array
            Parameters. Given as [a1,k1,a2,k2], the inputs into two von mises functions.
            
    Returns
    ---------------
    Array of firing rates given the var and par values
    '''
    f = pl.zeros(2*len(x))
    f[:len(x)] = ml.mises(x,var[0],par[0],par[1])*(1-ml.mises(x,var[1],par[2],par[3]))
    f[len(x):] = ml.mises(x,var[1],par[4],par[5])*(1-ml.mises(x,var[0],par[6],par[7]))
    return f

# set parameters
par  = [10,3,0.5,2,8,1.5,0.7,1]

# set bounds
bnds = ((-pi*1.5,pi*1.5),(-pi*1.5,pi*1.5))

# set starting conditions
var0 = [[-pi/2,-pi/2],[-pi/2,pi/2],[pi/2,pi/2],[pi/2,-pi/2]]

# number of neurons
N = 32

# basic x values for circular variables
A = pl.linspace(-pi,pi,N+1) #  in radiants
Ax = pl.linspace(-pi,3*pi,2*N+1)
Ax= Ax*180/pi # in degrees for plotting purposes

# number of contexts
nA = 16

# iterations 
iterations = 500

# to store estimations
a1_est = pl.zeros((nA,iterations))
a2_est = pl.zeros((nA,iterations))

# to store final biases
a1_bias = pl.zeros(nA)
a2_bias = pl.zeros(nA)

# stim 1
a1 = 0

# stim 2
A2 = pl.arange(-pi,pi,2*pi/nA)

# predefine holomaps
hmap_scatters    = hv.HoloMap(kdims=['context']) # the scatters
hmap_illusion1   = hv.HoloMap(kdims=['context']) # the illusion curve 1
hmap_illusionloc1= hv.HoloMap(kdims=['context']) # illusion dot 1
hmap_illusion2   = hv.HoloMap(kdims=['context']) # the illusion curve 2
hmap_illusionloc2= hv.HoloMap(kdims=['context']) # illusion dot 2
hmap_f           = hv.HoloMap(kdims=['context']) # population response
hmap_g           = hv.HoloMap(kdims=['context']) # stim 1 drive
hmap_h           = hv.HoloMap(kdims=['context']) # stim 2 modulation

# loop over stim 2's
for i,a2 in enumerate(A2):
    print 'Doing stimulus ' + str(i) + ' of ' + str(len(A2))
    # set variables
    var=[a1,a2]
    
    # get underlying firing rate
    f = fun(A,var=var,par=par)
    hmap_f[a2*180/pi] = hv.Curve(zip(Ax[:N],f[:N]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 1')
    
    # get second pop
#     g = fun(A,var=var,par=par1)
    hmap_g[a2*180/pi] = hv.Curve(zip(Ax[N:],f[N:]),vdims=[r_dim],kdims=[p_dim],label = 'Pop 2')
    
    # get modulation
#     h = fun(A,var=var,par=par2)
#     hmap_h[a2*180/pi] = hv.Curve(zip(Ax,h),vdims=[r_dim],kdims=[p_dim],label = 'Stim 2 modulation')
    
    # do iterations
    print 'iteration (max: ' + str(iterations) + '):'
    for it in range(iterations):
        print '\r'+str(it),
        # get noisy firing rate
        r = pl.poisson(f)
        
        # do the maximization
        L_max,var_max = ml.maxL(r,var0,bnds,fun,A,par,noise='poisson')
    
        # store variables
        a1_est[i,it] = var_max[0]
        a2_est[i,it] = var_max[1]
        
        # by convenction set sign properly (since can't be told by decoder)
        # which 
#         if pl.sign(a2_est[i,it]) != pl.sign(a2):
#             a2_est[i,it] *= -1
    
    # make sure everything is within -pi,pi
    a1_est[i,a1_est[i,:]<-pi]+=2*pi
    a1_est[i,a1_est[i,:]>pi] -=2*pi
    a2_est[i,a2_est[i,:]<-pi]+=2*pi
    a2_est[i,a2_est[i,:]>pi] -=2*pi
    
    # get bias
    a1_bias[i] = ml.circstats(a1_est[i,:])[0]
    a2_meanest = ml.circstats(a2_est[i,:])[0]
    
    # compare estimate for a2 with actual a2 to make circular
    a2_bias[i] = a2_meanest-a2#ml.circstats([a2_meanest,-a2])[0]
    a2_bias[a2_bias<-pi]+=2*pi
    a2_bias[a2_bias>pi] -=2*pi
    
    # fill in relevant holomaps
    hmap_scatters[a2*180/pi]    = hv.Scatter(zip(a2_est[i,:]*180/pi,a1_est[i,:]*180/pi),kdims=[a2_dim],vdims=[a1_dim]) # the scatters
    hmap_illusionloc1[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a1_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    hmap_illusionloc2[a2*180/pi]= hv.Scatter(zip([a2*180/pi],[a2_bias[i]*180/pi]),kdims=[a2_dim],vdims=[bias_dim])(style={'s':100}) # illusion dot 1
    
# fill in illusion curves
for i,a2 in enumerate(A2):
    hmap_illusion1[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a1_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 1
    hmap_illusion2[a2*180/pi]   = hv.Curve(zip(A2*180/pi,a2_bias*180/pi),kdims=[a2_dim],vdims=[bias_dim]) # the illusion curve 2
Doing stimulus 0 of 16
iteration (max: 500):
499 Doing stimulus 1 of 16
iteration (max: 500):
499 Doing stimulus 2 of 16
iteration (max: 500):
499 Doing stimulus 3 of 16
iteration (max: 500):
499 Doing stimulus 4 of 16
iteration (max: 500):
499 Doing stimulus 5 of 16
iteration (max: 500):
499 Doing stimulus 6 of 16
iteration (max: 500):
499 Doing stimulus 7 of 16
iteration (max: 500):
499 Doing stimulus 8 of 16
iteration (max: 500):
499 Doing stimulus 9 of 16
iteration (max: 500):
499 Doing stimulus 10 of 16
iteration (max: 500):
499 Doing stimulus 11 of 16
iteration (max: 500):
499 Doing stimulus 12 of 16
iteration (max: 500):
499 Doing stimulus 13 of 16
iteration (max: 500):
499 Doing stimulus 14 of 16
iteration (max: 500):
499 Doing stimulus 15 of 16
iteration (max: 500):
499
In [34]:
%%opts Curve {+axiswise}
# do holomap with scatter plots, showing them come closer as you move context
showfig  = hmap_illusion1*hmap_illusionloc1
showfig += hmap_illusion2*hmap_illusionloc2
showfig += hmap_f*hmap_g#*hmap_h
showfig += hmap_scatters
In [35]:
showfig.cols(2)
Out[35]:

Conclusion/discussion

  • Multiple variables make decoding/coding a lot harder due to ambiguities
  • Biases and decoding distributions can be very non-trivial
  • Could make predictions about bias distributions in humans

Thanks to...

  • Mark van Rossum, supervisor
  • Holoviews

Notes to self

  • multiple variables are okay, if other variables are for same preferred variable across population